# ICEMGenerator - 高效的ICEM CFD网格划分Replay Script生成工具
# 模块： Surface类，用于创建面，包括平面和曲面
# 作者： weixc
# 邮箱： 12028045 [at] zju.edu.cn
# 项目主页： https://gitee.com/weixinchong/icemgenerator

import numpy as np
from collections.abc import Iterable
from .point import Point
from .curve import Curve, Line, Ring, lines_from_points
from .param_check import check_axis

class SurfaceMethodError(Exception):
    '''面的构建方法错误'''
    def __init__(self, msg):
        super().__init__(msg)

class TooMuchEdgeCoexistError(Exception):
    '''客面与切割主面的重合线过多错误'''
    def __init__(self, ss, so, n_co_lines):
        ss_expr = str(ss)
        if ss.name != None:
            ss_expr += '(%s)' % ss.name
        so_expr = str(so)
        if so.name != None:
            so_expr += '(%s)' % so.name
        msg = '客面%s与主面%s重合的线为%d条，而切割过程中最多只允许出现一条重合的线' % (so_expr, ss_expr, n_co_lines)

class Surface():
    '''创建一个面'''

    curr_sId = 0  #类静态变量，用于确定最新创建面的编号
    curr_tmpId = 0  #类静态变量，用于确定切割平面时，产生的临时边缘线的编号
    members = []  #现有面的成员列表
    members_cut = []  #被切割的主面的成员列表

    def init_sId(sId: int):
        '''设置面的起始编号，应在开始创建点之前调用
        :param sId: 创建的第一个面的编号'''
        Surface.curr_sId = sId

    def create(curves: Iterable, name: str = None, method: str = 'round'):
        '''创建一个面（已存在则引用，不存在则创建）
        :param curves: 构成面的线，当面的构建方法`method`取如下不同值时，curves具有不同的含义:
            (1)round: 创建的面由封闭的线族curves围成，线无先后顺序要求，可用于创建平面、曲面；
            (2)sweep: curves内包括多条线，curves[0]为拉伸方向的轴线（可以是直线或弧线、圆环）后续的一条至多条线为基线，即创建以curves[1]~curves[n]为底，向curves[0]拉伸得到的曲面，可用于创建圆柱体、弯管等；
            (3)loft: curves内包括两条线，创建异面曲线curves[0]和curves[1]围成的曲面，线无先后顺序要求，可用于创建圆柱体、圆台等
        :param name: 面的名称，非必须
        :param method: 构建面的方法，即上面所述的round、sweep和loft
        :return (Surface): 创建的面'''
        # 执行参数检查
        if method not in ['round', 'sweep', 'loft']:
            raise SurfaceMethodError('面的构建方法%s不存在，method只能是`round`、`sweep`或`loft`' % method)
        if method == 'sweep' and len(curves) < 2:
            raise SurfaceMethodError('当采用sweep构建曲面时，在给定一条拉伸方向轴线的基础上，还应给出至少一条基线')
        elif method == 'loft' and len(curves) != 2:
            raise SurfaceMethodError('当采用loft构建曲面时，只能给定两条异面曲线')
        # 开始创建
        srf = Surface.find_member(curves)
        if srf == None:
            srf = Surface(curves, name, method)
        elif name != None:
            srf.setName(name)
        return srf

    def members_replay(cut_surfaces_priority: int = 0) -> str:
        '''生成所有面的Replay Script
        :param cut_surfaces_priority: 需要切割的面在创建时的优先顺序:
            (1)取0时，按原有创建面的顺序进行创建；
            (2)取1时，将这些面放在其他面之前进行创建；
            (3)取-1时，将这些面放在其他面之后进行创建
        :return: 生成的Replay Script'''
        if cut_surfaces_priority == 0:
            text = ''
            for s in Surface.members:
                text += s.replay() + '\n'
            return text
        else:
            text_merge = ''
            for s in Surface.members_cut:
                text_merge += s.replay() + '\n'
            text_no_merge = ''
            for s in Surface.members:
                if s not in Surface.members_cut:
                    text_no_merge += s.replay() + '\n'
            if cut_surfaces_priority == 1:
                return text_merge + text_no_merge
            elif cut_surfaces_priority == -1:
                return text_no_merge + text_merge

    def find_member(curves: Iterable):
        '''查找指定的面是否已经定义过
        :param curves: 要查找的面内的线
        :return (Surface | None): 找到已定义的面时，返回该面；否则为None'''
        for s in Surface.members:
            if len(curves) != len(s.curves):  #当两个面的线数量不一致时，两个面显然不是同一个面
                continue
            surface_is_same = True
            for c in curves:
                if c not in s.curves:
                    surface_is_same = False
                    break
            if surface_is_same:
                return s

    def __init__(self, curves, name, method):
        '''创建一个面'''
        self.sId = Surface.curr_sId
        if method == 'sweep':  #面构建方法为sweep时，曲线的顺序不能变
            self.curves = curves
        else:
            self.curves = sorted(curves, key = lambda c: c.cId)
        self.name = name
        self.method = method
        self.srf_cut = []  #被用于切割当前面中的面
        # 将各条线的点整合到面的点中
        self.points = []
        for c in self.curves:
            for p in c.points:
                if p not in self.points:
                    self.points.append(p)
        self.points = sorted(self.points, key = lambda p: p.pId)
        # 判断是否为平面，若是平面，则计算其法向量
        if method == 'round':
            if curves[0].type == 'ring':  #圆环围成的面即为平面
                self.is_plane = True
                self.direction = curves[0].direction
            else:  #直线围成的面
                self.is_plane = True  #先假定面为平面
                norm_first = None
                nc = len(curves)
                for i in range(nc - 1):
                    for j in range(i + 1, nc):
                        normv = np.cross(curves[i].direction, curves[j].direction)
                        if np.linalg.norm(normv) == 0:  #零向量不能用于判断方向，应跳过
                            continue
                        if norm_first is None:
                            norm_first = normv  #记录首个法向量
                        else:
                            if np.linalg.norm(np.cross(norm_first, normv)) > 0:  #根据法向量是否共线判断是否为平面
                                self.is_plane = False
                                break
                if self.is_plane:
                    self.direction = np.abs(norm_first) / np.linalg.norm(norm_first)  #对于平面，给定其法向量
        else:  #通过sweep、loft获得的面作为曲面处理
            self.is_plane = False
        Surface.members.append(self)  #将当前面增加到成员列表
        Surface.curr_sId += 1  #最新面编号增加1

    def __str__(self):
        return 'srf.%03d' % self.sId

    def __repr__(self):
        return self.__str__()

    def setName(self, name: str):
        '''更改面的名称
        :param name: 要设定的面名称'''
        self.name = name

    def replay(self) -> str:
        '''生成面的Replay Script
        :return: 生成的Replay Script'''
        nc = len(self.curves)
        if self.method == 'round':
            # 生成创建面的指令
            crv_expr = ' '.join([str(ci) for ci in self.curves])
            if nc >= 2 and nc <= 4:
                text = 'ic_surface 2-4crvs GEOM %s {0.01 {%s}}' % (self, crv_expr)
            else:
                text = 'ic_surface bsinterp GEOM %s {%s}' % (self, crv_expr)
            # 生成使用客面的边缘线切割主面的指令
            for obj, remain_cut_srf in self.srf_cut:
                crv_obj = obj.curves
                co_lines = find_co_lines(self, obj)  #重叠的线不能用于切割
                if len(co_lines) > 1:
                    raise TooMuchEdgeCoexistError(self, obj, len(co_lines))
                crv_expr = ''
                tmp_expr = ''
                i = 0
                for c in crv_obj:
                    if c in co_lines:
                        continue
                    if i == 0:
                        crv_expr += '%s' % c
                        tmp_expr += 'tmp%03d' % Surface.curr_tmpId
                    else:
                        crv_expr += ' %s' % c
                        tmp_expr += ' tmp%03d%d' % (Surface.curr_tmpId + (i - 1) // 10, (i - 1) % 10)
                    i += 1
                text += '\nic_geo_project_curve_to_surface {%s} %s %s GEOM 0 0' % (crv_expr, self, 'tmp%03d' % Surface.curr_tmpId)
                text += '\nic_geo_trim_surface %s {%s} 1' % (self, tmp_expr)
                if not remain_cut_srf:
                    text += '\nic_delete_geometry surface names %s.cut.0' % self  #删除切割后产生的多余部分
                Surface.curr_tmpId += 1  #切割完毕后更新临时边缘线编号
            return text
        elif self.method == 'sweep':
            script = ''
            for i in range(1, nc):
                if i >= 2:
                    script += '\n'
                script += 'ic_geo_cre_srf_crv_drv_srf GEOM %s %s %s 1' % (self, self.curves[i], self.curves[0])
            return script
        elif self.method == 'loft':
            return 'ic_geo_cre_srf_loft_crvs GEOM %s 0.1 {%s %s} 4 0 1' % (self, self.curves[0], self.curves[1])

    def cutBy(self, obj, remain_cut_srf: bool = False):
        '''使用客面切割当前面
        :param obj(Surface或Iterable<Surface>): 用于切割当前面的客面，可以是一个或多个
        :param remain_cut_srf: 是否保留切割后的客面'''
        if type(obj) == Surface:
            obj = (obj,)
        for oi in obj:
            self.srf_cut.append((oi, remain_cut_srf))
            Surface.members.remove(oi)  #切割后，客面应被移除
            if oi == self:  #如果客面跟当前面是同一个面，则当前面也应被移除
                Surface.members.remove(self)
        if self in Surface.members and self not in Surface.members_cut:
            Surface.members_cut.append(self)  #将当前面添加到需要切割操作的面的列表中

    def info(self, single: float = False) -> str:
        '''获得面的信息
        :param single: 是否输出单一信息，为True时仅输出面的信息，为False时还会输出每条线、每个点的信息
        :return 面的信息'''
        # 先输出面的信息
        text = 'Surface %03d' % self.sId
        if self.name != None:
            text += ' [%s]' % self.name
        if self.method == 'sweep':
            crv_expr = ' '.join([str(ci) for ci in self.curves[1:]])
            text += ': %s ~ (%s)' % (self.curves[0], crv_expr)
        else:
            crv_expr = ' '.join([str(ci) for ci in self.curves])
            text += ': (%s)' % crv_expr
        if not single:
            # 输出线的信息
            text += '\n[Curves Included]'
            for c in self.curves:
                text += '\n' + c.info(single = True)
            # 输出点的信息
            text += '\n[Points Included]'
            for p in self.points:
                text += '\n' + p.info()
        return text

    def getCurve(self, name: str) -> Curve | None:
        '''获得面中指定名称的线
        :param name: 线的名称
        :return: 找到指定名称的线时，返回该线；否则返回None'''
        for c in self.curves:
            if c.name == name:
                return c

    def getPoint(self, name: str) -> Point | None:
        '''获得面中指定名称的点
        :param name: 点的名称
        :return: 找到指定名称的点时，返回该点；否则返回None'''
        for p in self.points:
            if p.name == name:
                return p

def surface_from_points(points: Iterable) -> Surface:
    '''由一系列点绘制闭合的面
    :param points: 围成面的点
    :return: 绘制的面'''
    p_for_lines = points + [points[0]]
    lines = lines_from_points(p_for_lines)
    srf = Surface.create(lines)
    return srf

def create_rectangle(u1: float, u2: float, v1: float, v2: float, w: float = 0, axis: str = 'z', name: str = None) -> Surface:
    '''创建一个矩形
    :params u1,u2,v1,v2: 表明矩形由(u1,v1)、(u2,v2)围成
    :param w: 矩形的轴方向坐标
    :param axis: 矩形的轴向，可选值为x、y、z
    :param name: 面的名称，非必须
    :return: 创建的矩形平面'''
    check_axis(axis)
    # 创建点
    if axis == 'x':
        pnt_lb = Point.create(w, u1, v1, 'x->left_bottom')
        pnt_rb = Point.create(w, u2, v1, 'x->right_bottom')
        pnt_lt = Point.create(w, u1, v2, 'x->left_top')
        pnt_rt = Point.create(w, u2, v2, 'x->right_top')
    elif axis == 'y':
        pnt_lb = Point.create(u1, w, v1, 'y->left_bottom')
        pnt_rb = Point.create(u2, w, v1, 'y->right_bottom')
        pnt_lt = Point.create(u1, w, v2, 'y->left_top')
        pnt_rt = Point.create(u2, w, v2, 'y->right_top')
    elif axis == 'z':
        pnt_lb = Point.create(u1, v1, w, 'z->left_bottom')
        pnt_rb = Point.create(u2, v1, w, 'z->right_bottom')
        pnt_lt = Point.create(u1, v2, w, 'z->left_top')
        pnt_rt = Point.create(u2, v2, w, 'z->right_top')
    # 创建线
    crv_left = Line.create(pnt_lb, pnt_lt, 'left')
    crv_right = Line.create(pnt_rb, pnt_rt, 'right')
    crv_bottom = Line.create(pnt_lb, pnt_rb, 'bottom')
    crv_top = Line.create(pnt_lt, pnt_rt, 'top')
    # 创建面
    srf = Surface.create((crv_left, crv_right, crv_bottom, crv_top), name)
    return srf

def is_circle(srf: Surface) -> bool:
    '''判断一个面是否为圆面
    :param srf: 要判断的面
    :return: 是否为圆面'''
    return srf.curves[0].type == 'ring'

def get_uv(point: Point, axis: str = 'z') -> tuple:
    '''获得特定轴向的点的u、v的值
    :param point: 要读取的点
    :param axis: 点所在面的轴向，可选值为x、y、z
    :returns: u, v'''
    check_axis(axis)
    if axis == 'x':
        return point.X[1], point.X[2]
    elif axis == 'y':
        return point.X[0], point.X[2]
    elif axis == 'z':
        return point.X[0], point.X[1]

def get_uv_maxmin(srf: Surface, axis: str = 'z') -> tuple:
    '''获得一个面内各点u、v的最小、最大值
    :param srf: 要获取的面
    :param axis: 面的轴向，可选值为x、y、z
    :returns: u_min, u_max, v_min, v_max'''
    check_axis(axis)
    u_min = np.inf
    u_max = -np.inf
    v_min = np.inf
    v_max = -np.inf
    for p in srf.points:
        u, v = get_uv(p, axis)
        if u < u_min:
            u_min = u
        elif u > u_max:
            u_max = u
        if v < v_min:
            v_min = v
        elif v > v_max:
            v_max = v
    return u_min, u_max, v_min, v_max

def check_inclusion(srf1: Surface, srf2: Surface, axis: str = 'z') -> bool:
    '''判断客面是否被包含在主面内
    :param srf1: 主面
    :param srf2: 客面
    :param axis: 两个面的轴向，可选值为x、y、z
    :return: 客面是否在主面内'''
    check_axis(axis)
    if not is_circle(srf1) and not is_circle(srf2):  #主客面均为方形的情形（判断正确性不能保证）
        u1_min, u1_max, v1_min, v1_max = get_uv_maxmin(srf1, axis)
        u2_min, u2_max, v2_min, v2_max = get_uv_maxmin(srf2, axis)
        return u1_min <= u2_min and u1_max >= u2_max and v1_min <= v2_min and v1_max >= v2_max
    elif is_circle(srf1) and is_circle(srf2):  #主客面均为圆形时的情形
        u1, v1 = get_uv(srf1.points[0], axis)
        u2, v2 = get_uv(srf2.points[0], axis)
        l = np.sqrt((u1 - u2)**2 + (v1 - v2)**2)
        r1 = srf1.curves[0].r
        r2 = srf2.curves[0].r
        return l + r2 <= r1
    elif is_circle(srf1) and not is_circle(srf2):  #主面为圆形，而客面为方形的情形
        u1, v1 = get_uv(srf1.points[0], axis)
        r1 = srf1.curves[0].r
        for p in srf2.points:
            u2, v2 = get_uv(p, axis)
            l = np.sqrt((u1 - u2)**2 + (v1 - v2)**2)
            if l > r1:
                return False
        return True
    elif not is_circle(srf1) and is_circle(srf2):  #主面为方形，而客面为圆形的情形
        u2, v2 = get_uv(srf2.points[0], axis)
        r2 = srf2.curves[0].r
        for p in srf1.points:
            u1, v1 = get_uv(p, axis)
            l = np.sqrt((u1 - u2)**2 + (v1 - v2)**2)
            if l < r2:
                return False
        return True

def coordinates_all_same(values: Iterable, c1: int, c2: int) -> bool:
    '''检查所有点在c1、c2方向的坐标值是否均一致，以判断两条线是否重叠
    :param values: 要比较的坐标值的列表
    :params c1, c2: 方向的索引
    :return: 坐标值是否均一致'''
    v1 = values[0][c1]
    v2 = values[0][c2]
    for i in range(1, len(values)):
        if values[i][c1] != v1:
            return False
        if values[i][c2] != v2:
            return False
    return True

def find_co_lines(srf1: Surface, srf2: Surface) -> list:
    '''查找客面哪些线与主面重叠（重叠线的长度可以不一样）
    :param srf1: 主面
    :param srf2: 客面
    :return: 客面中与主面重叠线的列表'''
    co_lines = []
    lines_in_srf1 = list(filter(lambda c: c.type == 'line', srf1.curves))
    lines_in_srf2 = list(filter(lambda c: c.type == 'line', srf2.curves))
    for l1 in lines_in_srf1:
        for l2 in lines_in_srf2:
            if l2 is l1:  #两条线完全一致，则必然重叠
                co_lines.append(l2)
                continue
            p11 = l1.points[0].X
            p12 = l1.points[1].X
            p21 = l2.points[0].X
            p22 = l2.points[1].X
            if coordinates_all_same((p11, p12, p21, p22), 0, 1):  #点x,y坐标一致时，两条线重叠
                co_lines.append(l2)
            elif coordinates_all_same((p11, p12, p21, p22), 0, 2):  #点x,z坐标一致时，两条线重叠
                co_lines.append(l2)
            elif coordinates_all_same((p11, p12, p21, p22), 1, 2):  #点y,z坐标一致时，两条线重叠
                co_lines.append(l2)
    return co_lines

def create_circle(r: float, x: float = 0, y: float = 0, z: float = 0, axis: str = 'z', name: str = None) -> Surface:
    '''创建一个以axis为轴，半径为r的圆
    :param r: 圆的半径
    :param x, y, z: 圆心坐标
    :param axis: 圆的轴向，可选值为x、y、z
    :return: 创建的圆面'''
    check_axis(axis)
    # 创建线（圆环）
    edge = Ring.create(r, x, y, z, axis, 'ring')
    # 创建面
    srf = Surface.create((edge,), name)
    return srf